Center manifold reduction for large populations of globally coupled phase oscillators 
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A bifurcation theory for a system of globally coupled phase oscillators is developed based on the theory of 
rigged Hilbert spaces. It is shown that there exists a finite-dimensional center manifold on a space of generalized 
functions. The dynamics on the manifold is derived for any coupling functions. When the coupling function is 
sin 9, a bifurcation diagram conjectured by Kuramoto is rigorously obtained. When it is not sin 9, a new type of 
bifurcation phenomenon is found due to the discontinuity of the projection operator to the center subspace. 



Introduction. Collective synchronization phenomena are 
observed in a variety of areas, such as chemical reactions, en- 
gineering circuits, and biological populations Jl]]. In order 
to investigate such phenomena, a system of globally coupled 
phase oscillators of the following form is often used: 
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where 0i(f) denotes the phase of an z'-th oscillator, w, e R 
denotes its natural frequency drawn from some distribution 
function g(u>), K > is the coupling strength, and f(6) = 
Y^=-oo fi>e is a 27r-periodic function. When f(0) = sin 6, 
it is referred to as the Kuramoto model [2]. In this case, it is 
numerically observed that if K is sufficiently large, some or 
all of the oscillators tend to rotate at the same velocity on av- 
erage, which is referred to as synchronization lU [51] . In order 
to evaluate whether synchronization occurs, Kuramoto intro- 
duced the order parameter r(t)e ^^W, which is given by 
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When a synchronous state is formed, r(f) takes a positive 
value. Indeed, based on some formal calculations, Kuramoto 
assumed a bifurcation diagram of r(t): Suppose yV — > oo. If 
g(a>) is an even and unimodal function such that g"(Q) + 0, 
then the bifurcation diagram of r(t) is as in FigOJa). In 
other words, if the coupling strength K is smaller than K c := 
2/{7rg(Q)), then r(t) e is asymptotically stable. If K ex- 
ceeds K c , then a stable synchronous state emerges. Near the 
transition point K c , r is of order 0((K - K c ) l/2 ). See 0] for 
Kuramoto's discussion. 

In the last two decades, several studies have been performed 
in an attempt to confirm Kuramoto's assumption. Daido |4] 
calculated steady states of Eq. (Q~|) for any / using an argument 
similar to Kuramoto's. Although he obtained various bifurca- 
tion diagrams, the stability of solutions was not demonstrated. 
In order to investigate the stability of steady states, Strogatz 
and Mirollo and coworker performed a linearized anal- 
ysis. The linear operator T\, which is obtained by linearizing 
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FIG. 1. Bifurcation diagrams of the order parameter for (a) f(9) = 
sin 9 and (b) f(S) = sin 9 + h sin 29. The solid lines denote stable 
solutions, and the dotted lines denote unstable solutions. 



the Kuramoto model, has a continuous spectrum on the imag- 
inary axis. Nevertheless, they found that the steady states can 
be asymptotically stable because of the existence of resonance 
poles on the left half plane [8]. Since the results of Strogatz 
and Mirollo and coworker are based on the linearized anal- 
ysis, the effects of nonlinear terms are neglected. However, 
investigating nonlinear bifurcations is more difficult because 
T\ has a continuous spectrum on the imaginary axis, that is, a 
center manifold in the usual sense is of infinite dimension. In 
order to avoid this difficulty, Crawford and Davies [9] added 
a noise of strength D > to the Kuramoto model. The con- 
tinuous spectrum then moves to the left side by D, and thus 
the usual center manifold reduction is applicable. However, 
their method is not valid when D — 0. An eigenfunction of T\ 
associated with a center subspace diverges as D — > because 
an eigenvalue on the imaginary axis is embedded in the con- 
tinuous spectrum as D — > 0. Recently, Ott and Antonsen ifioll 
found a special solution of the Kuramoto model, which allows 
the dimension of the system to be reduced. Their method is 
applicable only for f(ff) = sin 9 because their method relies on 
a certain symmetry of the system 111 ill . Furthermore, the re- 
duced system is still of infinite dimension, except for the case 
in which g(a>) is a rational function. Thus, a unified bifurca- 
tion theory for globally coupled phase oscillators is required. 

In the present letter, a correct center manifold reduction 
is proposed by means of the theory of generalized func- 
tions, which is applicable for any coupling function /. It 
is shown that there exists a finite-dimensional center mani- 
fold on a space of generalized functions, despite the fact that 
the continuous spectrum lies on the imaginary axis. This 
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will be demonstrated for two cases, (I) f(6) = sin 6 and (II) 
f(6) = sin 8 + h sin 26, h e R, and two distribution functions 
g(oS), (a) a Gaussian distribution and (b) a rational function 
(e.g., Lorentzian distribution g(oj) = 1/(tt(1 + w 2 ))). For (I), 
we rigorously prove Kuramoto's conjecture, and for (II), a dif- 
ferent bifurcation diagram will be obtained, as was predicted 
by Daido Q4f] . The different bifurcation structure is shown to 
be caused by the discontinuity of the projection to the gener- 
alized center subspace. All omitted proofs are given in 111211 . 
Settings. The continuous model of Eq. (Q~|) is given by 

dpr/dt + d(p t v)/d6 = with v := to + ^ZS-<x,/aWe" V=Tw . 
where 77/ is defined to be 

JrJo 

and p, = p t {0,a>) is a probability measure on [0,27r) pa- 
rameterized by t, to e R. In particular, 771 is a continuous 
version of Kuramoto's order parameter. Setting Zj(t,to) : = 

J^e ^ i6 p,{6, to)d6 yields 



By the canonical inclusion, Eq. (0) is rewritten as an evolution 
equation on the dual space X' as 



j t I Zj) = T* \Zj}+ ^PljK J] MP \Z,)-\ Z H ), (6) 



where T* is a dual operator of Tj defined through (<p \ Tjp) = 

Here, the strategy for the bifurcation theory of globally cou- 
pled phase oscillators is to use the space of generalized func- 
tions X' rather than a space of usual functions. The reason for 
this is explained intuitively as follows. If we use the space 
L 2 (R, g(to)dto) to investigate the dynamics, then the behavior 
of p, itself will be obtained. However, it is neutrally stable 
because of the conservation law C p,(8, a>)d6 = 1 . What we 
would like to know is the dynamics of the moments of p t , in 
particular, the order parameter. This suggests that we should 
use a different topology for the stability of p,. (Note that the 
definition of stability depends on definition of the topology.) 
For the purpose of the present study, p, is said to be convergent 
to p as t — > 00 if and only if 
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(t,to)g(to)dto = (Zj,P ) = (P ,Zj), 
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where ( • , • ) is the inner product on the weighted Lebesgue 
space L 2 (R, g(ui)du>), and Po(oj) = 1. The linear operator 
Tj is known to have a continuous spectrum on the imagi- 
nary axis. Furthermore, there exists a positive constant IQr, 
such that if K^r < K, Tj has eigenvalues on the right half 
plane (such that the de-synchronous state is unstable), while 
if < K < K { J\ Tj has no eigenvalues. For example, if / 
is an odd function and if g is even and unimodal, then IQr — 
-\m(fj)l(n\fj\ 2 g(§)). In the present letter, for simplicity, we 
assume that K c := inf, K ( c f> = (for f(ff) = smd + h sin 20, 
which is true if and only if h < 1). When < K < K c , Tj 
has no eigenvalues, and thus the dynamics of the linearized 
system dZj/dt - TjZj is quite nontrivial. In II 1211 . the spec- 
tral theory on rigged Hilbert spaces is developed to reveal the 
dynamics of the linearized system. 

A rigged Hilbert space consists of three spaces X c 
L 2 (R, g((L>)da>) c X': a space X of test functions, a Hilbert 
space L 2 (R, g{u)doj), and the dual space X' of X (a space 
of continuous anti-linear functionals on X, each element of 
which is referred to as a generalized function). We use Dirac's 
notation, where for fi e X' and cp e X, fj.(<p) is denoted by 
(0 1 p). For a e C, we have a(<p \ p) = (5(p \p) = (<p \ ap). The 
canonical inclusion i : X — > X' is defined as follows. For 
ifr G X, we denote i(tfr) by \if/), which is defined to be 



for any /' e Z and (p € X. The topology induced by this conver- 
gence is referred to as the weak topology. By the completion 
of L 2 (R, g(aj)dcj) with respect to the weak topology, we ob- 
tain a space of generalized functions X' . On the space X', a 
function Z\(t, at) converges as t — > 00 if and only if {(p\Z\) 
converges for any <p e X. Since the order parameter is written 
as 771(f) = (Po I Zj), this topology is suitable for the purpose of 
the present study. For this topology, it turns out that p, is not 
neutrally stable. 

A suitable choice of X depends on g{to). When g(a>) is a 
Gaussian distribution, let Exp + (/3) be the set of holomorphic 
functions defined near the upper half plane such that \(p(z)\e~® z ' 1 
is finite. Set X = Exp + := U/?>o Exp + (/8). We can introduce 
a suitable topology on Exp + so that the dual space Exp^ be- 
comes a complete metric space, which allows the existence 
of a center manifold on Exp^ to be proven. When g(a>) is a 
rational function, X := H+ is a space of bounded holomor- 
phic functions on the real axis and the upper half plane. In 
this case, we can show that i(H+) c H' + is a finite-dimensional 
vector space, which implies that Eq. © is essentially a finite- 
dimensional system. This is why in flCl [1311 . the system is 
reduced to a finite-dimension system when g(oj) is the (sum 
of the) Lorentzian distribution. 

Let e Ti ' be the semigroup generated by T\ . In lfl2Tl . the spec- 
tral decomposition of e T[t is obtained by means of the rigged 
Hilbert space. Define resonance poles Aq, A\, ■ ■ ■ of T\ to be 
roots of the equation 
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/((A)(0) = <0 1 iff) := (<f>, .A) = <Ka>Ma))g(a>)da). (5) 



on the imaginary axis and the left half plane. Roughly speak- 
ing, a resonance pole is a continuation of an eigenvalue to the 
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second Riemann sheet of the resolvent (A - T{)' A HQ. In 
the following, for simplicity, we assume that all A n 's are single 
roots of Eq. (Q. Define a functional fij e X' to be 



R Aj - V-Tfc 
+2*0(- V=U;)*(- V=T^), (Re(^) < 0), 
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The yu/s are called the generalized eigenfunctions associated 
with the resonance poles due to the equality T*\fij) = Aj\ fij). 
Then, we can prove that the semigroup is expressed as 
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for any ifr e X, which gives the spectral decomposition of e T,t 
on the dual space X', where D n are constants defined by 
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(In the previous paragraph, X was chosen so that the right- 
hand side of Eq.© converges.) Equation © completely de- 
termines the dynamics of the order parameter for the lin- 
earized system. If < K < K c , then all resonance poles 
A„ lie on the left half plane. As a result, 771(f) = (e T, '<f),Po) = 
(cp I (e T,t ) x Po) decays to zero exponentially as f — > 00, which 
proves the asymptotic stability of the de-synchronous state. 

Center manifold reduction. When K = K c , there exist 
resonance poles on the imaginary axis. The generalized cen- 
ter subspace E c c X' is defined as a space spanned by gener- 
alized eigenfunctions associated with resonance poles on the 
imaginary axis, say Ao, ■ ■ ■ , A M - Equation (© suggests that the 
projection IT C to E c is given by 



M K 

IW> = V— (^l^-l/Un). 
n=0 AUn 
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In general, He : X' — > X' is continuous only on a subspace of 
X' because the topology on X' is too weak. When X = Exp + , 
it is proven in IU2I1 that II C is continuous only on /(Exp + (0)). 
For a solution of Eq. d3j, Zi,Z2, - are included in Exp + (0), 
although Z_i, Z_2, - are not. Because of the discontinuity of 
II C , an interesting bifurcation occurs when f(8) + sin#. In 
what follows, assume that g is an even and unimodal function. 
In this case, on the imaginary axis, there is only one resonance 
pole /lo = at K = K c . Hence, E c = span{//o) is of one di- 
mension, where // ( > is the generalized eigenfunction associated 
with Aq = 0. Next, let us derive the dynamics on a center man- 
ifold. The derivation is performed in the same way for both 
(a) a Gaussian distribution and (b) rational functions. 

(I) First, we consider f(ff) - sin 8. In this case, equations 
for Z\ , Z2, - do not depend on Z \ , Z_2, • • • • Thus, II C acts 



continuously on solutions of Eq. d6). Set s — K - K c . Then, 
Eq. (O for j — 1 is rewritten as 

^|Zi> = rfoIZO + |(P |Zi>|Po> - |<P |Zi>l^>, (10) 

where Tio is defined by replacing K with K c in the definition 
of T\ . In order to obtain the dynamics on the center manifold, 
using the spectral decomposition, we set 

|Z 1 ) = ^a| A i > + |y 1 ), a(t) = ^-(Z l \ f x ), (11) 

along the direct sum E c © E^. The purpose here is to de- 
rive the dynamics of a. Since | Y\) and \Zq) are included in 
the stable subspace, the center manifold theorem [ 12] reveals 
that on the center manifold, | Y\), \ Z2) ~ 0(a 2 ). In particular, 
the last term <P I Z t >| Z 2 ) of Eq. QU]) is of order <9(a 3 ). Ap- 
ply the projection IT to both sides of Eq. ( TTOb - Noting that 
<Po |Zi)II c |Z2) is also of order 0(a 3 ) because the projection 
is continuous on | Zq), we obtain the dynamics on the center 
manifold as 



a 
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dt DoK c 
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If g"(Q) < 0, this equation has a fixed point of order 0((K - 
K c ) 1/2 ) when s = K - K c > 0. It is easy to verify that the 
order parameter 771(f) is given by 771(f) = a + h.o.t. Thus, the 
dynamics of the order parameter is also given by Eq. (TT2V 
Since Dq > 0, when g is even and unimodal, a = (de- 
synchronous state) is unstable, and the nontrivial fixed point 
(synchronous state) is asymptotically stable when s = K - 
K c > 0, which confirms Kuramoto's diagram. 

(II) Assume that f(ff) = sin6» + h sin 26 with h € R. Then, 
Eq. (|6]l for j — 1 is given by 

j t \Z 1 ) = T* \Z 1 ) + ^(P \Z l )\P ) 

- I «P 1 Zi>| Z 2 > + h(P 1 Z 2 )| Z 3 > - h(P 1 Z 2 )| Z_!» . (13) 

In this case, Z_i, on which Il c is discontinuous, appears. As 
before, \Za) satisfies \Zq) ~ 0(a 2 ) and ndZz) ~ 0(a 2 ) since 
II C acts continuously on IZ2). This implies that the term 
(Pq |Zi)|Zz) in Eq. ( TT3l yields a cubic nonlinearity as case 
(I). On the other hand, since II C is discontinuous on Z_i, we 
can show that II t | Z_i) ~ 0(1) even if | Z_i) ~ O(a). As a re- 
sult, the last term (Pq \ Zq)\ Z_i) in Eq. (TT3l yields a quadratic 
nonlinearity. Indeed, we can verify that 

rT | Z_!> = ^£^ e -V=Targ W + 0(a) (14) 

Do 

Applying the projection Il c to both sides of Eq. ( fT3l . we ob- 
tain the dynamics on the center manifold as 

± = _2L L - &Lae-+***) + h.o.t., (15) 
dt D K C \ 2(1 -h) ) K ' 
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FIG. 2. Numerical results for f{9) = sin 9 + h sin 29. Black 
dots denote the order parameter calculated from Eq. ([T] for N = 
8,000,g(w) = e-" 2 ' 2 / V2tt using the method shown in 1 14]. Since r 
is unstable when K — K c < 0, it is difficult to obtain small r. The solid 
lines are interpolations of black dots using quintic polynomials. The 
dotted lines denote the analytical results obtained by Eq. d!6t . 



where C := p.v. Lg'(tL>)/coda> is a negative constant, which 
proves that there exists a fixed point that is expressed as 



l^il ~ M = ^^T (K " + 0((K ~ ^ 

K,.Cn 



Therefore, for h < 0, a stable branch emerges when K c < K, 
and for < h < 1 , an unstable branch emerges when K < K c 
(Figffib)). 

Discussion. Equation (O shows that the dynamics of 
Z\ , 7a, • ■ ■ is independent of Z \ , Z 2, - if and only if f(ff) = 
sin 9. In other words, Eq. © splits into two systems: a system 
of {Z\,Zq, • ■•} and a system of {Z-\,Z-z, ■ ■ ■}. Since the pro- 
jection Il c is continuous on a solution of the former system, 
we can show the existence of a smooth center manifold. Note 
that Eq. ([1]) is invariant under the rotation on a circle. As a re- 
sult, the dynamics on the center manifold is also invariant un- 
der the rotation ffiHe ^~^a. If a center manifold is smooth, 
then the dynamics on this manifold with the rotation symme- 
try must be of the form a = aF(\a\ 2 ). Thus, a cubic nonlinear- 
ity is dominant, and a pitchfork bifurcation generally occurs, 
as shown in Eq. ( fTSl i. On the other hand, if f(0) + sin#, 
then the equations of Z\,Zq,, ■ ■■ depend on Z \,Z i, • •• , on 
which Il c is not continuous. In such a case, the center mani- 
fold is not smooth, and quadratic nonlinearity may appear, as 
described above. In this manner, different bifurcations occur 
when f{9) + sin 9. Although the diagram shown in FigJUb) 
looks like a transcritical bifurcation, Eq. (Q3) is different from 
the normal form of a transcritical bifurcation. Because of the 
factor e~^ mg ^ caused by the discontinuity of Il c , Eq. ( fT3] l 
remains invariant under the rotation despite the existence of a 
quadratic nonlinearity. The discontinuity induces a new type 
of bifurcation including <?~ V-Targ(<r)_ 

A center manifold reduction for globally coupled phase os- 
cillators was also developed by Crawford and Davies H with 
a noise of strength D > 0. Although they also expected a 



diagram such as shown as FigJTfb) when D = 0, the factor 
e~ V-Targ(a) was not Stained. Since the eigenfunction diverges 
as D — > 0, expressions of the dynamics on the center manifold 
were not shown explicitly. In the present letter, we have shown 
that the eigenfunction fio exists on a space of generalized func- 
tions, which provides a correct center manifold reduction. The 
diagram shown in FigQlb) was also obtained by Daido [4] by 
means of a self-consistent analysis. Unfortunately, his results 
were not correct because he performed inappropriate termwi se 
integrations of certain infinite series. According to his results, 
the order parameter is given as ( 1 - 2h) ■ const., which suggests 
that some degeneracy occurs when h = 1/2. However, the 
numerical results given in Fig|2]show that the critical expo- 
nent of the order parameter changes only when h — 1 , which 
agrees with the results of the present study (IT6V Ott and An- 
tonsen [10] found an inertia manifold given by Z„ = (Z\) n 
when f(9) = sin 9. The center manifold of the present study 
is a finite-dimensional submanifold of the inertia manifold, 
which provides a further reduction of the results of Ott and 
Antonsen. The key strategy of the present theory is to use 
spaces of generalized functions and the weak topology. The 
weak topology is suitable for investigating the dynamics of 
moments of probability density functions. Since the strategy 
is independent of the details of the models, this strategy will 
be extended to various types of large populations of coupled 
systems and evolution equations of density functions, such as 
the Vlasov equation. 
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